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T— I Abstract 

o 

^-H Approximate Bayesian Computation has been successfully used in population ge- 

netics models to bypass the calculation of the likelihood. These algorithms provide 
an accurate estimator by comparing the observed dataset to a sample of datasets 
^ simulated from the model. Although parallelization is easily achieved, computation 

times for assuring a suitable approximation quality of the posterior distribution are 
still long. To alleviate this issue, we propose a sequential algorithm adapted from 
Del Moral ct al. (2012) which runs twice as fast as traditional ABC algorithms. Its 
parameters are calibrated to minimize the number of simulations from the model. 
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1 Introduction 



We work in the parametric Bayesian setting and are interested in the posterior distribution 
of e £ e C W^. Let TT{e) denote the prior distribution and /(x|0) the hkehhood. The 
probabihty density function of the target posterior is given as 

7r(0|Xobs) oc 7r(0)/(xobs|^) , 

where Xobs G P is the observed dataset, not necessarily an independent and identically 
distributed (iid) sample. For complex models, the computation of values of the likelihood 
function /(xobs| ■ ) can be expensive or plain impossible which renders sampling from 
■ l^obs) extremely difficult. Approximate Bayesian Computation (ABC) is a recent 
technique that only requires sampling from /( • |0) (Marin et al., 2012). This method- 
ology comes from the population genetics community. We refer the reader to Marin 
et al. (2012) and Beaumont (2010) which provides an excellent survey of ABC methods 
and its uses in biology. For more practical considerations, one can also consult Csillery 
et al. (2010). The likelihood- free rejection sampler of Pritchard et al. (1999) works as 
follows: 

for i = 1 to do 
repeat 

Generate 6' from the prior distribution '7r( • ) 
Generate z from / ( • |0') 

untild(5(z),S(xobs)) <e 
end for 

We simulate first a parameter value 0' from the prior, then a dataset z from the corre- 
sponding model. We subsequently compare the simulated dataset to the observed one 
Xobs via some summary statistics S :T) ^ S and a distance function d : S x S ^ M^. If 
z is close enough to Xobs according to a tolerance level e, we keep 6' . The likelihood- free 
rejection algorithm samples from 

7r,(0|xobs) oc /. ^ ^ 7r(0)/(z|0)dz, (1) 

J|zGC|rf(5(z),S{xob,))<ej 

the marginal distribution of 7re(0,z|xobs)- The idea behind ABC is that the summary 
statistics coupled with a small tolerance should provide a good approximation of the 
posterior distribution 7re(0|xobs) ~ 7r(0|xobs)- Marjoram et al. (2003) have introduced a 
Monte Carlo Markov Chain (MCMC) algorithm targeted at (1) which does not require 
any calculation of the likelihood. 



2 



Rejection sampling and MCMC methods can perform poorly if the tolerance level e is 
small. In practice, the tolerance level e used in the rejection sampling algorithm is not fixed 
in advance, but corresponds to a quantile of the distances between the observed dataset 
and some simulated ones. Various sequential Monte Carlo algorithms (Doucet et al., 
2001; Del Moral, 2004; Liu, 2008) have been constructed as an alternative to these two 
methods, see Sisson et al. (2007), Sisson et al. (2009), Beaumont et al. (2009), Drovandi and 
Pettitt (2011) and Del Moral et al. (2012). These algorithms start from a large tolerance 
level Eq, and at each iteration the tolerance levels decreases, et < £t-i- The simulation 
problem becomes therefore more and more difficult, whereas the proposal distribution for 
the parameters 6 becomes more and more accurate. 

The algorithm of Beaumont et al. (2009) corrects the bias of the one of Sisson et al. 
(2007) (Sisson et al., 2009) and is a particular Population Monte Carlo scheme (Cappe 
et al., 2004). It requires fixing a sequence of decreasing tolerance levels Eq > ei > . . . > et 
which is not very realistic for practical problems. In contrast, the proposals of Del Moral 
et al. (2012) and Drovandi and Pettitt (2011) are adapted likelihood-free versions of the 
Sequential Monte Carlo sampler (Del Moral et al., 2006) and include a self-calibration 
mechanism for the sequence of decreasing tolerance levels. 

We consider here the specific case where all the others calculations are negligible com- 
pared to a random generation from the model. This is typically true for complex models, 
e.g. for complex scenarios in population genetics. We adapt the likelihood-free SMC sam- 
pler of Del Moral et al. (2012) such that the number of simulated values from the model 
is as small as possible. 

We present a naive likelihood-free sequential algorithm and we recall the proposals of 
Del Moral et al. (2012) and Drovandi and Pettitt (2011) in Section 2. Then we introduce 
our self-calibrated strategy in Section 3. Finally, we illustrate its numerical behavior on 
simulated data and a challenging real-data example from population genetics. 



2 Sequential likelihood-free schemes 

As already explained, the likelihood-free rejection scheme consists in generating a sample 
|(0j,Xj), 1 < i < from 7r(0)/(x|0) and in keeping the 0i's such that 

(i^5(xj) , S'(xobs)^ < We obtain the collection {^j}j^jf^^j^y where J(e,A^) = |j G 

{l,...,iV} : Xj G {z e V\d{S{z),S{xohs)) < e}}- If e is too small, the set J{e,N) is 
almost empty. On the contrary, if e is too large then 7r£(0|xobs) can be very far from 
the posterior distribution. In practice, a quantile a of the distances d( 5(xj), 5(xobs) ) is 
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chosen which corresponds to setting e such that J(e,A^) contains [aA'^J distinct indices 
(where [aj denotes the integer part of a). With such a strategy, the 0j's are always 
sampled from the prior distribution. This can be very inefficient, especially when vr(0) 
is not informative. Sequential schemes aim at constructing proposals for the parameters 
which iteratively gain relevance. We start by introducing a naive approach based on the 
SMC sampler of Del Moral et al. (2006). 

From now on, we assume a prior '7r(0) that is uniform over the compact set 6prior- 
In an initial step, all sequential strategies employ the likelihood-free rejection sampler to 
create a sample from the distribution vrej, (0|xobs) • In the further steps, a population of 

N particles := |(0*,x*), 1 < i < drawn from the distribution 7rej(0|xobs) is given 
at iteration t + 1. The naive strategy at iteration t + 1 is based on a fixed sequence of 
quantiles {at}(>o ("^o associated to Eq) and a given target value for e : 

1. Sort the particles |(0*,x*),l < i < with respect to their distances from Xobs, 
i.e. 

d(^Si4),S{^obs)) < d(5(x*),5(x„fe,)) < ••• < d(^S{^%),Si^obs)) 

2. Calibrate £t+i- 

Set Et+i such that et+i = ci(^5(x[^^^_^jyj), S'(xobs)) 

3. Compute the new particle system: 

(a) Create a new set of particles V* in the following way: 

- Repeat [l/a^+ij times the first [at+iA^J particles of 

- Draw the N — / at-\-i\at-^-iN \ last particles of V* randomly among 

{(0*,x*),l<i< L«mA^j} 

(b) Move each particle (0*,x*) of V* according to a MCMC kernel Kt with sta- 
tionary distribution iTet+i (^j z|xobs) to get a new set of particles 7^*+^ 

4. l( Et+i < e, return "P^+i and et+i as output of the algorithm. Otherwise, set t = 
and return to step 1. 

Moving according to Kf. 

We propose to use a Metropolis-Hastings kernel in step 4. For instance, we can com- 
pute the empirical variance of the {0*}'s in the set Vf For each particle (0*,x*) of 
V*, we generate some 6 according to the Gaussian distribution J\f[6*,2T^^ and some x 
according to the likelihood /(x|0) . Since the prior is uniform, the acceptance rate for such 
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a Gaussian proposal is binary, see Del Moral et al. (2012) for more details. If G Sprior 
and d(^S(x),5(xobs)) < et+i, we set (0-+\x*+^) = (^,x). Otherwise, we do not move 
and set x*+i) = (0*, x*) . 

No Rao-Blackwellization. 

We might replace Step 3. (a) by the following alternative: Given "P*, the system V* 
consists of iid particles generated from a uniform distribution over the [afA^J first par- 
ticles of 7^*. If iV — [[l/tttJatA'^J = 0, Step 3. (a) corresponds to a conditional expectation 
under this alternative: we integrate over the uniform draws from the [oiA^J first particles 
of P*. 

Quality of the output. 

A sequential schemes returns weighted particles (xi,(x)i), . . . , (x^r, a;Ar) associated to 
a target distribution Q such that 

N 

1=1 

for any Q— measurable set A. One way to measure the quality of a weighted sample is 
the Effective Sample Size (ESS) criterion, a quality indicator inspired from importance 
sampling ideas (cf. (Liu, 2008, chapter 2)). It is given by 

^ -1 

ESs((xi,^i),...(x^,6.^)) = (^^^2^ . (2) 

i=l 

In many situations, the particle system is obtained by a resampling step followed by a 
systematic assignment of the weight 1/N to each particle. Such a system can be composed 
of a small number of distinct particles repeated many times. Then we can not directly 
apply the formula (2). Indeed, ESS would equal A^. In this case, we calculate the ESS 
by considering only the distinct particles with their aggregated weight: if one particle 
X* is repeated r times with the weight 1/A^, its weight after agregation is r/N. If there 
are repetitions, another indicator of the quality of the weighted sample is the number 
of distinct particles. If the particles are repeated the same number of times, it is equal 
to the ESS. This indicator is an upper bound for the ESS. The difference between these 
two measures is even greater when the entropy of the weights is strong. If no movement 
is accepted at the step 3 of the above-described naive strategy, the final particle system 
consists of aoai . . . utN distinct particles with the same weight. Each one is repeated 
many times in the final system. In other words, this is equivalent to running the likelihood- 
free rejection sampler with a quantile a equal to a^ai . . .ax- We thus obtain an ESS that 
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corresponds to the number of distinct particles, 

ESSt = aooi • • • otN, 

which can indicate very poor quaUty. If now we assume that at the iteration t a proportion 
pt of particles is accepted in the step 3 (b) of the naive strategy, the number of distinct 
particles is 

ao{ai + pi){a2 + P2) ■ ■ ■ {ar + Pt)N. (3) 

To explain this value of ESS, we assume that at iteration t the system contains Nt 
distinct particles and that at = with an integer divisor it of Nt. If each particle 
(0, x) proposed during step 3 is accepted in the system with probability pt, then the 
system V* consists of it copies of the [a* A^J first particles of "P*. Hence the mean number 
of distinct particles in the system 7^*+^ is 

et 

atNt + ^ atptNt = atNt + ItatPtNt = {at + pt)Nt. 

i=l 

The likelihood-free SMC proposals of Del Moral et al. (2012) and Drovandi 
and Pettitt (2011). 

The proposal of Del Moral et al. (2012) is an adaptive scheme that calibrates the 
threshold £t by controlling the ESS during the iterations. The threshold £t is chosen such 
that 

ESSt+i = aESSt, 

where a e]0,1[. In practice, they consider the number of particles in the iteration t for 
which the distance to the observation is lower than £t+i- This quantity is treated as the 
ESS and repetitions are not taken into account. The parameter a is a quality index. 
If a ~ 1 then the scheme will move slowly towards the target ABC and the successive 
approximations are good. Conversely, if a is close to 0, the stabilization at the desired level 
is very fast but the quality of the successive approximations is poor. In their numerical 
experiments, they keep a between 0.9 and 0.99. This scheme is not optimized to limit the 
number of simulations according to the model. Typically, for the same parameter value 
they propose to draw M samples from the augmented ABC target 

M M 
k=l k=l 

Finally, we note that they perform a resampling step before the kernel shifting step only 
if the ESS falls below a certain threshold. 
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The proposal of Drovandi and Pettitt (2011) applies the naive strategy by fixing the 
sequence of quantiles at to a unique value 1 — a. Only the [aN\ particles with the largest 
distance from the observation are modified. This scheme is also clearly not calibrated to 
minimize the number of simulations according to the model. 



3 Our self-calibrated algorithm 

The main idea of our method is to choose the minimal et+i for which (a^ + pt) > 1. In the 
above algorithm, we then have to permute step 3 (in which et+i is fixed) and step 4 (in 
which we can compute pt). Unfortunately the MCMC kernel Kt of step 4 depends heavily 
on et+i, hence this permutation is not trivial. Similar to the naive approach, we assume 

that at iteration t + 1 we are given a population of N particles := | (0*, x*) ,1 < i < 

drawn from the distribution ir^^ (0|xobs)- At iteration t + 1, our self-calibrated likelihood- 
free SMC sampler works as follows: 

1. Sort the particles |(0-,x*), 1 < i < with respect to their distances to Xq^js) i-S- 
di = d(s(x*),5(xobs)) <d2 = d(5(x*),5(xobs)) <---<dN = d(s(x*v),5(xobs)) 

2. Calibration of et+v 

Compute cr|, the variance of {s-^Si < i < iv} in the particle system V . 

a' ^ 0.0 
repeat 

a' ^a' + 0.01 

e{a') ^ d^^'ivj 

for i = 1 to [a'iVj do 

if particle {6i,^ij does not exist then 

Generate Oi from the normal distribution J\f {6\,2a1^ and Xj from /(x|0j) 

end if 
end for 

Set N' {a'^ as the number of particles (Oi^^i) among the [a'A^J first ones that 

satisfy 6i £ Opnor and d(^S{xi),S{xohs)^ < e'(a') 

p'{a') ^ N'{ay[a'N\ 
until a' + p'{a') > 0.9 
at+i ^ a', et+i ^ e' {a') and pt+i ^ p'(a') 
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3. Computation of the new particle system: 

for i = to [of+iA^J do 

if 9i G Gprior and d(^S{Sci),S(xohs)^ < £t+i then 

else 

end if 
end for 

for i = [at+iN\ + 1 to A/" do 

Generate an integer number / uniformly between 1 and [at+iATj ^ 
Generate 6 from the normal distribution M{6\, ^crf) and x from /(x|0) 

if e Gprior and (i(^5(x) , S'(xobs)) < £t+i then 

else 

end if 
end for 

4. If pt+i < 0.1, return 7?*+^ and et+i as output of the algorithm. Otherwise, set 
t = t + l and return to step 1. 



Adaptive scheme and final tolerance level. 

Clearly, p' {a') increases with a' since the conditions 6i E Q and d^S'(xi) , S'(xobs)^ < 

e'(a') become less restrictive and the Metropolis-Hastings kernel moves a larger proportion 
of particles. Hence a' + p' [a') increases with a'. 

Two important questions need closer examination: 



1. What is the behavior of the threshold of acceptance e at each iteration when N 
tends to infinity? 

2. What can we say about the final level of acceptance in our scheme? 

When N increases, the calibrated tolerance level et does not converge to 0, but to 
a positive quantity that depends only on the prior it, the likelihood / and the transi- 
tion mechanisms Kt during the Metropolis-Hastings step. At each stage of the itera- 
tive algorithm, the triplet (a, p, e) takes values on a one-dimensional manifold. If we 
parametrize this manifold by e, the calibration scheme chooses the tolerance level e such 
that a(^s) + pi^s) = 1. 
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Consider for instance the {t + l)th iterative stage of the algorithm. When N ^ oo, 
the empirical distribution of the particle system converges to the joint distribution 
TTej (0|xobs)/(x|0). Hence asymptotically a(e) is the cumulative distribution function of 

(i( S'(x) , S'(xobs) ) when (0,x) is drawn from VTe^ (0|xobs)/(x|0) • In other words, 



Q(e) = / / 7ret(0|xobs)/(x|0)l r , , ^ dxd0, for any e > 



The rejection step of the {t + l)th iterative stage keeps a proportion a(e) of the 0's. 
Hence, the empirical distribution of those 0's converges to 

7r£(0|xobs) = a(e) W vret(0)/(x|0)l j- ^ dx, for any e > 0. 

j c((^5(x),s(xob,)j<£ ^ 

Asymptotically, the acceptance rate p[e) for the Metropolis-Hastings step associated to 
the Markov kernel i^t+i is 

p{e)= fll Tie{0\^o^.s)Kt+i{e*\e)f{^*\e*)i. ^dxWde. 

d(^5(x*),5(xobs)j<£ 

Therefore the calibrated tolerance level £t+i is asymptotically the solution of 

dxd0 



//vr,,(0|xobs)/(x|0)lr y 
•' •' <^ d( 5(x),s(xob=) )<e > 




+ / // ^,(0|xobs)i^m(0ie)/(x*|r)i . ^ dx*drd0 = 1 

d( s(x*),5(xobs) )<£ 



In turn, calibrated parameters of the other iterative stages are consistent estimators of pos- 
itive quantities that depend only on the prior, likelihood and the kernel inside Metropolis- 
Hastings. A full and formal proof should show the convergence of these empirical quanti- 
ties. Calibrated parameters converge when one reaches the asymptotic regime with large 
N (see numerical experiments on figures 1 and 2). The moment of reaching the stopping 
criterion [p < 0.1) does not vary significantly across independent runs of our scheme. This 
convergence is illustrated on the toy example of Sisson et al. (2007) in Section 4.1. 



9 



For the analysis of the quahty of the final threshold e when the number of particles A'^ 
increases to infinity, we assume that the final level of acceptance of our ABC-SMC scheme 
does not vary significantly from one instance to another. The ABC approximation given by 
the output level et (where T is the stopping time) may not be low enough to assure a good 
approximation of the posterior. Actually, e is often too large and we have to add an extra 
rejection step on the output of the sequential algorithm to decrease e. Adding an extra 
rejection step is faster than violating the stopping criterion and increasing the number of 
iterations. In addition, the ESS of the resulting simulated sample of parameter-dataset 
pairs depends linearly on (before and after the extra rejection step). 

Remember that the ABC approximation vTe is reconstructed from the final sample 
via classical methods (kernel density estimate, local linear regression,. . . ). We can infer 
its characteristics such as the mean, the median or quantiles with classical sample-based 
estimators. The reconstruction error depends heavily on the ESS of the simulated sample. 
Increasing increases the quality and accuracy of the reconstruction. When a good 
ABC approximation is desired, we can distribute the computational cost by running the 
instances of our ABC-SMC scheme in parallel and merging the outputs. The level of 
acceptance of the resulting sample is determined as the greatest level among the distributed 
samples. 

High Performance Computing. 

A major advantage of the classical ABC algorithm is its aptitude for parallelization on 
multicore architectures and computing clusters. In our algorithm, we could conduct the 
simulations according to the probabilistic model in parallel. However this demands for a 
shared memory architecture to run the sequential algorithm since each iteration depends 
on the previously simulated sample. If the simulation of a dataset from the model and 
the computation of the summary statistics are slow, we observe a good parallel program 
performance on a multicore computer using the OpenMP API (see http://openmp.org): 
The algorithm runs almost eight times faster when using eight cores than when using a 
single core. 

The situation is different if one wants to parallelize the program on several nodes of 
a cluster. The overhead due to communication and synchronization between the nodes 
might be substantial, with the trickiest part being the calibration of the tolerance level at 
each step. In this case, we recommend independent runs of the algorithm on each node of 
the cluster, providing independent outputs we can easily merge. In particular, we assume 
that we are able to construct random number generators on each node of the cluster which 
are mutually independent. 
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Efficient learning. 



The aim of iterative algorithms is to learn progressively how to draw whithout 
introducing bias into the computation of the posterior. Indeed, simulating 6 from the prior 
distribution can lead to simulated data sets far from the observation. In our proposal, we 
learn the current ABC approximation of the posterior from the previously simulated pairs 
(0, x) and we draw pairs (6, x) from this current approximation in order to gradually 
decrease the tolerance level e. Clearly, applying the learning scheme is inefficient if the 
posterior distribution is almost equal to the prior distribution (see the discussion on the 
toy example of Sisson et al. (2007) in section 4.1). By the way, we can interpret the 
Metropolis-Hastings step as drawing 0' from a kernel density estimate of the current ABC 
approximation of the posterior. Actually, we draw the parameter 0' of the proposed couple 
(0',x') from the Gaussian mixture 



[atN\ 



where ht = y^2cj^ is the smoothing parameter. This mixture is a smooth density estimate 
of 1^*"^, 1 < i < [a^A^J I whose empirical distribution is a proxy for vTe^. 

Comparison vi^ith the classical ABC algorithm. Our ABC-SMC scheme returns 
a final particles system |(0f,x^),l < i < a| with ESSy and final acceptance level 
£t- We consider this final output as an iid simulated sample of size ESSt- The time 
complexity of our algorithm in the setup of population genetics is determined by the total 
number of datasets simulated from the model. We compare the tolerance level reached by 
our algorithm with the one reached by the classical acceptation-rejection method while 
demanding the same computational effort and the same quality of the output in terms of 
ESS. We denote by Nt the total number of simulations according the model required for 
our scheme. The classical ABC algorithm with the same time complexity draws Nt pairs 
(0,x) from 7r(0)/(x|0). The output of this classical algorithm should be of size ESSt 
to be comparable to the output of our proposal. Hence one should compare £t with the 
quantile £e,cc-re]{ESST/NT) of order ESSt/Nt of D = d{S (x) , S {xobs)) when (0,x) is 
drawn from the joint 7r(0)f{x\6). To quantify the efficiency of our scheme compared to 
the acceptation-rejection algorithm, we compute the gain factor. We compute the number 
of particles A'acc-rej needed when using the ABC acceptation-rejection sampler to provide 
an approximation with level £t and ESSt particles accepted. The gain factor of our 
ABC-SMC scheme in terms of the number of simulations according to the model is given 

by 

rr=^. (4) 
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In terms of simulation according to the model, our scheme outperforms the standard 
acceptance-reject algorithm when this gain factor is greater than 1. 

Convergence properties. 

If the decreasing sequence of tolerance levels {et}j>Q is fixed, the successive target 
distributions {'^st}t>o random. If moreover the Markov kernels {-RTtj^^Q are 

deterministic, the convergence results on the SMC algorithm of Del Moral et al. (2006) 
can be applied directly. 

If we consider the case where the Markov kernels depend on the past particles, we can 
provide convergence results using the recent works of Del Moral et al. (2011). We can 
also use theorems on triangular arrays of random variables given in Douc and Moulines 
(2008). Therefore, denote by ■ = {e%i,x%,^,l < i < the system of N particles 
at the end of the t-th iteration of the algorithm. If we introduce the sub-cj-field J'^^ = 
^{^N'h 1 ^ ^ ^ the particles j, I < i < A^j are independent given Tj^^. Under 
standard regularity conditions, we can prove the following law of large numbers and CLT: 



If f |(^(^) Ivr^j (d^) < +00, the convergence in probability 

1 ^ 

i=l 

holds with TTej ((^) := f ipi^^'^n^^ (d^) as tends to infinity. 
If / y?^ (.^)7rej (d.^) < cxD, the convergence in distribution 

^ ^AA(0,5?(v.) 

holds with S^[(p) := / (/?^(^)7r£j (d^) — 7r£^((/?)^ as N tends to infinity. 



4 Numerical experiments 
4.1 A toy example 

The following toy example, first introduced in Sisson et al. (2007) and widely used in 
the literature, consists of a simple mixture of two Gaussian distributions with different 
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variances and unknown mean: 



-10,10 ' 

11 (5) 
/(x|0) = -0(x; e, 1) + -</)(x; e, 1/100) , 

where ^a,fe] denotes the uniform distribution on the interval [a, 6] (a < 6) and 0(-;/u,cr^) 
the probabihty density function of the one-dimensional Gaussian distribution with mean 
H and variance o"^. As explained in Beaumont et al. (2009) and Del Moral et al. (2012), 
the posterior distribution associated to the observation Xobs = is such that 



7r(6>|xobs) oc |0(6>;O,1) +(/.(6>; 0,1/100) I Ir 

J ■^-10<6»<10> 

This example is simple enough to compare the iterative methodologies. 



(6) 



Table 1 provides a comparison of our ABC-SMC scheme to the ones proposed by 
Del Moral et al. (2012) and Drovandi and Pettitt (2011) in terms of the total number of 
simulations according to the model necessary to provide an approximation with the same 
acceptance level e = 0.09. The schemes proposed in Del Moral et al. (2012) and Drovandi 
and Pettitt (2011) are not designed to optimize the number of simulations according to 
the likelihood. Clearly, our algorithm performs best. 



Table 1: Comparison of computational cost of the different ABC-SMC schemes. 





Cost 


ESS 


Final acceptance level e 


Del Moral et al. (2012) 


46 X 10^ 


29250 


0.09 


Drovandi and Pettitt (2011) 


109 X 10^ 


32000 


0.09 


Our ABC-SMC scheme 


23 X 10^ 


33285 


0.09 



We have applied our proposal with 10^ particles at each iteration. Figures 1, 2 corre- 
spond to the 1st and 8th iterative stage and illustrate the convergence of the calibrated 
tolerance levels e when the simulated sample size increases. We also observed that the 
variance of e is proportional to 1/A^. 

We mentioned earlier that sequential schemes are not recommended when the prior 
provides sufficient information about the parameter. Figure 3 (two top curves) shows 
that when the prior density support is close around 0, the stopping criterion is met in 
the early iterations and the gain factor is not favorable to the use of sequential schemes. 
The four curves in this figure show the evolution of the gain factor introduced in (4) 
over the iterations. The vertical dotted line (blue) indicates the stop time. In each of 
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Figure 1: Left: confidence intervals of probability 0.95 for the calibrated tolerance level 
e at the end of the first iterative stage when N increases. Right: inverse of e's variance 
with respect to N. Variances and confidence intervals were computed on 250 independent 
replications. 




Figure 2: Left: confidence intervals of probability 0.95 for the calibrated tolerance level 
e at the end of the 8th iterative stage when increases. Right: inverse of e's variance 
with respect to A^. Variances and confidence intervals were computed on 250 independent 
replications. 
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the four studies, we initialize our scheme with a particle system obtained by acceptance- 
rejection with a quantile level uq. Thus, we conclude that our ABC-SMC scheme reduces 
significantly the number of simulations according to the model when the prior provides 
little information (see the bottom curves in figure 3). 

When the prior density support is very large, the sequential scheme needs many itera- 
tions to detect the relevant region for the particle simulation. This can be very inefficient 
and may have a strong negative impact on the gain factor. The choice of the quantile level 
ao is important in this situation. A good choice of ao avoids the region where learning 
is useless. The first ABC approximation vTep (0|xobg) assigns a very low mass to those 
parts of the prior density support. We illustrate this choice of ao in figure 3 (bottom 
curves). With prior density support [ — 100,100] and ao = 1/4, the maximum of the 
gain factor is 14. This choice of ao is obtained using a system of particles drawn from 
the joint distribution 7r(0)/(0|x). The level ao is chosen such that the variance of Oi 
selected after one stage of acceptance-rejection is significantly reduced. More specifically, 
the choice ao = 1/4 reduces the standard deviation of by a factor 4. 

4.2 A population genetics example 

We now consider an experiment that relates more directly to the genesis of ABC, namely 
population genetics. We propose to perform a Bayesian inference of the parameters of an 
evolutionary scenario. This scenario retraces the invasion of bees in Europe. The scenario 
is presented in Figure 4. It consists of five inter-populational events which occurred at 
instants ti, . . . ,1^ (from the oldest to the newest). It is composed of two types of events: 
three divergences parameterized by their dates ti,t2, ti and two admixtures parameterized 
by their dates ^3 , and rates ri and r2 . 

This scenario involves six populations that have diverged at different moments, as 
described in Figure 4. It includes two unobserved populations. We use different colors to 
distinguish the evolution of the different populations in the Figure 4. 

The datasets for this study are composed of different samples of individuals collected 
from four populations observed until present and denoted by Popl, . . . , PopA in Figure 4. 
Each sample is of the order of fifty. Our datasets consist of individual genotyping at eight 
independant microsatellite loci. To perform the ABC analysis, we consider 30 summary 
statistics selected from the statistics available in the DIYABC software of Cornuet et al. 
(2008). 

The GSM model is considered in the study. The mutation rate is given per unit of 
time and per individual. For the i**^ loci of the study, the total mutation Hi^cLO rate is 
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divided into two parts as follows: 



fJ'i,GLO = fJ'i,sm- 

is the rate that corresponds to an increase or a decrease of the DNA sequence by a 
length G X m, where G is a realisation of a geometric random variable with parameter 
0.42 and m is a pattern length. The rate /ij,sNi is associated to the addition or removal of 
one nucleotide. 

Priors. 

Denote by Nei, . . . , Ncq the effective population sizes associated to the six popula- 
tions involved in the scenario. This size varies with the population. Recall that inter- 
populational events in the scenario occurred at ti, . . . ,tQ. Thus, by regrouping the ad- 
mixture rates ri,r2 and the mutation rates /ij, /Uj^sNi we obtain a vector of parameters of 
dimension 15. 

The priors on the parameter vector are given as follows 

• We suppose that 

Nei ~ Gamma[o.i;5ooooo] {Ne/2, 2) , 

where Gamma^Q.^ (p, s) is the gamma distribution with position and shape parame- 
ters p, s truncated to [a;b]. We have E [A^Cj I A^e] = Ne, and A^e is an hyper-parameter 
with uniform law ^5;ioo,ooo]- 

• the priors associated to the different event dates are given as ~ ^o.Ol;800]l U~*5 ~ 

^0;50,000]i *3 - ^5 ~ ^[0;50,000]) ^2 - max (is, U) ~ ^100;500,000] ) ^1 - ^2 ~ ^100;2000,000] • 

• The priors on the mutations parameters are given as 

/i, ~M Gamma (7Z/2,2), 

where JI = 10~^ and E [/ij] = ]I. We also have 

IJ'i,SNi ~iid Gamma ^^^_g,^^_3j {jIsni/'^,'^), 

where Jlgj^j = 10"*^ and E[ni^sNi] = J^sNi- 

• The priors on admixture rates ri and r2 are uniform on [O.Ol; 0.999] . 

Figure 5 presents the results for the estimates of the parameters 9s defined for each 
population as 9 = 4A^e/i. One can see that the replication power of the introduced 
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algorithm is significative. Figure 6 shows the posterior approximations for the different 
admixture rates. FinaUy, the gain factor for our proposal compared to the standard ABC 
acceptation-rejection agorithm is given in Figure 7. The dashed vertical line indicates the 
stopping time of the iterative algorithm, i.e., the first T for which pT < 0.1. We stop 
the algorithm after 14 iterations. We only need half the number of simulations from the 
model to obtain results similar to the standard ABC acceptation-rejection algorithm. 



5 Discussion 



We proposed a self-calibrated ABC-SMC algorithm adapted from Del Moral et al. (2012) 
such that the number of simulations from the model is minimized. On a very challenging 
population genetics example, we shown that our proposal runs twice as fast as traditional 
ABC algorithms. 

Our ABC-SMC scheme stops when the acceptance rate of Hastings-Metropolis falls 
below a fixed threshold. It seems possible to estimate the gain factor rt and stop when 
it falls below 1. Typically, for each iteration, we must be able to estimate the probability 
that the distance to the observation is smaller than et 



where the distribution of x is the marginal of the joint distribution (0,x) ~ 7r(0)£(x|0). 
This is equivalent to estimating the number of simulations required for the ABC acceptance- 
rejection sampler to accept ESS* particles. We have 



^d(5(x),S(xobs)) < et) =P(d(5(x),5(xobs)) < St d(5(x), 5(xobs)) < St-i 

p((i(5(x),5(xobs)) <ei|f^(5(x),5(xobs)) <eo) 
p(d(5(x),5(xobs)) <eo) . 
Thus, we can envision to estimate the probability of the rare event 

|d(5(x),5(xobs)) <et^ 

by the product of the levels of quantiles at calibrated by our scheme. At iteration t, x is 
distributed according to the marginal of 7r(0)/(x|0) given (i^5(x) , 5(xobs)^ < ^t-i- It 
would be very interesting to analyze the theoretical properties of this criterion. 
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Figure 3: Study of the efficiency of our proposal based on the information provided by 
the prior. Top, left: the prior density support is [—0.1; 0.1] and ao = 1/10. Top, right: the 
prior density support is is [—1; 1] and ao = 1/10. Bottom, left: the prior density support 
is [—10; 10] and ao = 1/2. Bottom, right: the prior density support is [—100; 100] and 
ao = 1/4. 
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Figure 7: Time factor for the population genetics example 
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